#include "config.h"

#include <math.h>

#include <gsl/gsl_rng.h>

#include <clx/libECM.h>

#include "libms2sim.h"

typedef struct {
	gsl_rng	*rng;
	double	*probs;
	unsigned int	size;
	double	skew;
	double	bottom;
} popularity_t;

static double *
build_probabilities(unsigned int size, double skew)
{
	double	*probs;
	double	sum = 0;
	unsigned int	i;

	probs = (double *)ecm_malloc(sizeof(double) * size);
	if (probs == NULL)
		return NULL;
	for (i = 1; i <= size; i++) {
		probs[i - 1] = 1.0 / pow(i, skew);
		sum += probs[i - 1];
	}

	probs[0] /= sum;
	for (i = 2; i < size; i++)
		probs[i - 1] = probs[i - 2] + probs[i - 1] / sum;

	return probs;
}

msl_popularity_t
msl_popularity_create(unsigned int size, double skew, unsigned long seed)
{
	gsl_rng	*rng;
	popularity_t	*popularity;

	rng = gsl_rng_alloc(gsl_rng_default);
	if (rng == NULL)
		return NULL;
	popularity = (popularity_t *)ecm_malloc(sizeof(popularity_t));
	if (popularity == NULL) {
		gsl_rng_free(rng);
		return NULL;
	}
	if (seed > 0)
		gsl_rng_set(rng, seed);
	popularity->rng = rng;
	popularity->size = size;
	popularity->skew = skew;
	if (size > 1) {
		popularity->probs = build_probabilities(size, skew);
		if (popularity->probs == NULL) {
			gsl_rng_free(rng);
			ecm_free(popularity);
			return NULL;
		}
	}
	else {
		popularity->probs = NULL;
	}
	return popularity;
}

unsigned int
msl_popularity_get(msl_popularity_t msl_popularity)
{
	popularity_t	*popularity = (popularity_t *)msl_popularity;
	double	prob;
	unsigned int	i;

	if (popularity == NULL)
		return 0;

	prob = gsl_rng_uniform(popularity->rng);
	for (i = 1; i < popularity->size; i++) {
		if (prob < popularity->probs[i - 1])
			break;
	}

	return i;
}

void
msl_popularity_free(msl_popularity_t msl_popularity)
{
	popularity_t	*popularity = (popularity_t *)msl_popularity;

	if (popularity == NULL)
		return;

	gsl_rng_free(popularity->rng);
	ecm_free(popularity->probs);
	ecm_free(popularity);
}
